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Abstract. We determine the complete asymptotic behaviour of the work distribution in driven stochastic 
systems described by Langevin equations. Special emphasis is put on the calculation of the pre-exponential 
factor which makes the result free of adjustable parameters. The method is applied to various examples 
and excellent agreement with numerical simulations is demonstrated. For the special case of parabolic 
potentials with time-dependent frequencies, we derive a universal functional form for the asymptotic work 
distribution. 
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1 Introduction 

With the discovery of work [T] and fluctuation [2113] the- 
orems in stochastic thermodynamics (for recent reviews 
see [H[S]), the traditional emphasis of statistical mechanics 
on averages was extended to include also large deviation 
properties. Indeed, averages like the one appearing in the 
Jarzynski equality 



,-0F _ 



-pw\ 



(1.1) 



are dominated by unlikely realization of the random vari- 
able (here the work W), and detailed information about 
the tail of the corresponding probability distribution is 
necessary to obtain an accurate result (see e.g. [§]). By 
definition, rare realizations are difficult to get, and conse- 
quently, numerically and even more experimentally gener- 
ated histograms seldom reach far enough into the asymp- 
totic regime. It is therefore desirable to have some addi- 
tional and independent information about the asymptotic 
behaviour of the relevant probability distributions. 

In the present paper we use the so-called method of 
optimal fluctuation to determine the asymptotics of work 
distributions in driven Langevin system. Special emphasis 
is put on the calculation of the pre-exponential factor that 
makes any fitting between histogram and asymptotics su- 
perfluous. By considering various examples, we show that 
the results for averages like (jl.l|) improve significantly if 
the pre-exponential factor is included. We use a novel 
method [7] to determine this pre-factor which builds on 
the spectral ^-function of Sturm-Liouville operators. The 
method is very efficient and straightforward in its numer- 
ical implementation. It also allows to handle the case of 
zero modes which is relevant in the present situation. For 



harmonic potentials with time-dependent frequency, we 
are able to derive the general form of the asymptotics of 
the work distribution analytically. 

The paper is organized as follows. Section [5] gives the 
basic equations and fixes the notation. In section [3] we re- 
call the basic steps in the determination of functional de- 
terminants from spectral (^-functions and adapt the pro- 
cedure to the present situation. Section 0] discusses two 
examples, one that can be solved analytically and merely 
serves as a test of the method, and one for which the 
analysis has to be completed numerically. In section [5] we 
elucidate the particularly interesting case of a harmonic 
oscillator with time-dependent frequency. Here, substan- 
tial analytic progress is possible. Finally, section [6] con- 
tains some conclusions. Some more formal aspects of the 
analysis are relegated to the appendices A-C. 



2 Basic equations 

We consider a driven stochastic system in the time interval 
< t < T described by an overdamped Langevin equation 
of the form 



x = -V'(x,t) + 



(2.1) 



The degrees of freedom are denoted by x, the time-depen- 
dent potential V gives rise to a deterministic drift, and 
is a Gaussian white noise source obeying (£(£)) = 
and (£(£)£(£')) = S(t—1f). Derivatives with respect to x are 
denoted by a prime, those with respect to t by a dot. The 
system is coupled to a heat bath with inverse temperature 
p. The initial state x(t = 0) =: xq of the system is sampled 
from the equilibrium distribution at t = 



* Dedicated to Werner Ebeling on the occasion of his 75th 
birthday. 



Po(x ) = — exp(- pV (x )) 
^0 



(2.2) 
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with Vo(x) := V(x, t = 0) the initial potential and 

Z = [ dxexp(-0Vo(x)) (2-3) 



the corresponding partition function. 

During the process, the system is externally driven and 
the potential changes from Vo(x) to Vr(x) := V(x,t = T) 
according to a given protocol. The work performed by the 
external driving depends on the particular trajectory x(-) 
the system follows and is given by [5] 



W[x(-)} = / dtV(x(t),t) 
Jo 



(2.4) 



Due to the random nature of x(-), also W is a random 
quantity. According to the general rule of transformation 
of probability, its pdf is given by 



P (W)= J e -PV M f dxT 



x(T)=x T 

x J Dx(-) p[x(-)] 6{W - W[x{-)}) . 

x(0)=x o 

The probability measure in trajectory space is [5] 



(2.5) 



p[ aJ (.)]=^ex P (-|jf at (x + V\x,t)) 2 ) , (2.6) 

where for mid-point discretization in the functional inte- 
gral we have 



TV 



< x P( J j dt V"(x(t),tfj 



(2.7) 



Using the Fourier representation of the S- function in (|2.5p . 
we find 



P(W) =jvj^jdx 



dq 

4tt//3 



T 

x(T)=x T 



Vx(-) e-W'Ud (2.8) 



x(0)=x o 



with the action 



s\,v). q ] = i ;,(,-„ H j dt[±(x+v>) 2 +jV 



- l j W . (2.9) 



The asymptotic behaviour of P(W) may now be deter- 
mined by utilizing the contraction principle of large devi- 
ation theory [TIT]. Roughly speaking, this principle stipu- 
lates that the probability of an unlikely event is given by 
the probability of its most probable cause [TTUT2"] . In the 
present context this means that whereas typical values of 
W are brought about by a variety of different trajectories 
x(-), the rare values from the tails of P(W) are predom- 
inantly realized by one particular path x(-) maximizing 



P[x(-)} := po(xo)p[x(-)] under the constraint W = W[x(-)]. 
A convenient way to implement this idea in the present 
context, is to evaluate the integrals in (j2.8[) by the saddle- 
point method. This is formally equivalent to considering 
the weak noise limit j3 — > oo. 

Let us therefore study expression (|2.8p in the vicinity 
of a particular trajectory x(-) and a particular value q of 
q. We put x(t) = x{t) + y(t) and q = q + r and expand 
up to second order in y(-) and r. After several partial 
integrations, we find 



S[x(-), q] = S + Sy la + Squad + 



where 
S 

"Sim 



■S[x(-),q] 
1 



(2.10) 



(2.H) 



(V ' - x )y + (V£ + x T )yi 



dt (x + V' -V'V" -iqV')y 

o 



ir(W-j dtV 



(2.12) 



quad 



(VoVo - 2)0)2/0 + (VtVt + yr)yi 



L dt <-d? 



+ 2ir / dt V'y 



V" 2 + V'V"' ~{l~iq)V" )y 
(2-13) 



Here, the notation V := V(x(t),t) and similarly for the 
derivatives of V has been used. 

The most probable trajectory a;(.) realizing a given 
value W of the work is specified by the requirement that 
5nn has to vanish for any choice of y(-) and r. We hence 
get the Euler-Lagrange equation (ELE) 

x + (1 - iq)V' - V'V" = (2.14) 

together with the boundary conditions 

x - Vq = 0, xt ~\~ Vj^ — . (2.15) 

From the term proportional to r we find back the con- 
straint 



W 



dt V 



(2.16) 



Note that for given W the solution for x(-) is usually 
unique and includes the optimal choice of its initial and 
final point. 

Neglecting contributions stemming from Squad, we ar- 
rive at the estimate 



P(W) ~ exp(-/3S*) 



(2.17) 



giving the leading exponential term for the asymptotic 
behaviour of P(W). It is solely determined by the optimal 
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trajectory itself: Using the properties (|2.14l) . (|2.15[) and 
(|2.16l) of x(-) in (|2.9[) . one can show 



W 1 - - 1 
S = - Y + -(V T + Vo)- ~{x T V^ + x V{) 



1 



+ - / dtV \V' - xV") + 



1 



iq 



dtxV . (2.18) 



Assuming that all eigenvalues are strictly positive, as nec- 
essary for x(-) being a minimum of < S[a;( > )], the c„ integrals 
may be performed and we find 

/¥ 1 f dr ( (3r 2 ^d 2 \ ,„ „„, 



In many cases an improved estimate including the pre- 
exponential factor is desireable. To obtain it, also the 
neighbourhood of the optimal trajectory has to be taken 
into account. This is possible by retaining S'quad in the 
exponent and by performing the Gaussian integrals over 
yo,yr,y(-) and r. 

In order to calculate 



vTI 

Integrating finally over r, we are left with 
V2 1 



I 



(2.27) 



Using 



\ [ \„ = detA and X) i = {V^-^V) - ( 2 -28) 

"71 



y{T)=y T 



dyo / dy T 



dr 

4tt//3 

y(o)=yo 



Vy(-) e -/ 3S Wb(-).'i ( (2. 19 ) we finally get 

Z 



we determine the eigenvalues A„ and normalized eigen- 
functions <p n (i) of the operator A defined by (cf. (|2.13[1 ) 



P(W) 



detA^^A-^V) 



==(1+0(1/0)). (2.29) 



.4 



+ {V"f + V'V" - (1 - iq)V" (2.20) 

dt z 

together with the Robin-type boundary conditions 

V r / Vn(0)-^„(0) = 0, V%<p n (T) + tp n (T) = 0. (2.21) 
Next we expand 

y(t) =J2cn<Pn{t) (2-22) 



and replace the integrations over yo,yr and y(-) by inte- 
grations over the expansion parameters c n according to 



dc„ 



=p (- § E A " c » - '-r E <=»*•) ■ < 2 - 23 > 



Here, we have introduced the notation 



cL := 



dt ¥>«(*) V"(f) 



(2.24) 



o 



and J is a factor stemming from the Jacobian of the trans- 
formation of integration variables. With this transforma- 
tion being linear, the Jacobian is a constant; with trans- 
formations between the eigenfunction systems of different 
operators being orthogonal, this constant cannot depend 
on the special form of V. In appendix |21 we show by com- 
parison with an exactly solvable case that 



This is the main result of this section. The same expression 
was obtained in |13j from a discretization of the functional 
integral. 



3 Calculating determinants from spectral 
^-functions 

The determination of the extremal action S occurring in 
(|2T2^)) requires the solution of the ELE ([2"TT4")l . Although 
this can be done analytically only in a few exceptional 
cases, its numerical solution poses in general no difficulty. 
Similarly, the calculation of (F'|^4 _1 |V'') is rather straight- 
forward. One solves the ordinary differential equation 

A1>(t) = V'(t) (3.1) 
with boundary conditions (|2.2ip and uses 

(V'lA-^V)^ fdtmvXt). (3.2) 
Jo 

The determination of det A is somewhat more involved. 
We will use a method introduced recently [7] , see also [13] , 
building on the spectral C-fvmction of Sturm-Liouville op- 
erators. The essence of the method is contained in the 
relation 

det 4 F(0) 



detA Ie{ F rof (0) ' 



(3.3) 



where the eigenvalues A„ of A are given by the zeros of 
-FX A), and similarly F re f (A) = determines the eigenvalues 
of the reference operator A lc f . Moreover, one has to ensure 
F(A)/F re f(A) — > 1 for |A| — > oo. For operators of the type 
considered here, 



1- IE. 

J V f> 



(2.25) 



.4 



d^ 
~dt 2 



■9(t) 



(3.4) 
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with homogeneous Robin boundary conditions 

a ip(t = 0) + b (p(t = 0) = , (3.5) 
c <p(t = T) + d<p(t = T) = , (3.6) 

a convenient choice for F(X) is [7] 

F(X) = c X x(T) + d±x(T) , (3.7) 

where x\ (t) is the solution of the initial value problem 

A X x(t) = \xx(t), Xa(0) = -b, xa(0) = a . (3.8) 

Note that x\(t) is defined for general A. The initial condi- 
tion in (|3.8[) is chosen such that xa (t) satisfies the bound- 
ary condition (|3.5p at t = for all values of A. Only if 
A coincides with one of the eigenvalues of A, A = A n , 
the boundary condition at t — T is satisfied as well. In 
this case xa is proportional to the eigenfunction ip n cor- 
responding to A„. The roots of the equation F(X) = are 
therefore indeed the eigenvalues of A. 

In appendix [B] we show that for the reference operator 



A ref = - 



dt 2 



with boundary conditions (|3.5p . (|3.6p one finds 

detA rc f _ 2 
Frof(0) ~ ~bd ' 



(3.9) 



(3.10) 



Hence, combining (|3.3[) and (|3.10p . the determination of 
det A boils down to the solution of the initial value prob- 
lem (|3.8p . Numerically, this is again straightforward. 



4 Examples 

4.1 The sliding parabola 

The simplest example for the class of problems considered 
is provided by a Brownian particle dragged by a harmonic 
potential moving with constant speed |T5TT6 llT7l[T8] 



V(x,t) 



(x - tf 



(4.1) 



In this case, the full distribution P(W) is known analyti- 
cally [HE]: 



(4.2) 



(4.3) 

This example hence merely serves as a test of our method. 
To apply (|2.29p . we first note that 




N = e T ' 2 and Z Q = JtL 



(4.4) 



Moreover, the ELE (|2.14p is linear and can be solved an- 
alytically: 



x(t; W) = -{2t + e- f - e l - T ) - W 



(2 - e- 1 - e*- 1 ) 



2(T + e- T - 1) 
Using this result in (|2.18p . we find 

- (W~(T + e- T -l)) 2 



4(T + e- T - 1) 



(4.5) 



The operator A defined in (|2.20p is, for the potential (I4.ip . 
given by 

A = -^ + l (4.6) 

with boundary conditions 

<p n (0) - p„(0) = 0, <p n (T) + cp n (T) = . (4.7) 

In order to determine F(0), we have to solve the initial 
value problem 

-Xo + Xo(*) = 0, X o(0) = 1, Xo(0) = 1 . (4.8) 

The solution is xo(t) = e* implying (cf. ([37?)) F(0) = 2e T . 
Using p. 31) and (|3.10p we hence find 



det A = 4e T 



(4.9) 



Finally, V' = -1, and in order to calculate (V'\A- 1 \V'), 
we have to solve 

m - 4>(t) = i , 

^(0)-^(0)=0, V(r)+^(T) = (4.10) 
which gives 

1 



m= +e -2). 

Plugging this into p.2p yields 

(V'lA-^V') = T + e- T -1 



(4.11) 



(4.12) 



Combining (j2~2^1) . (|47o]l , gU), and (j432l . we find 

back (|4.2I) . Since P(IU) is Gaussian, the quadratic expan- 
sion around the saddle-point already reproduces the com- 
plete distribution, i.e. there are no higher order terms in 



4.2 The evolving double-well 

As a more involved example, we discuss the time-depen- 
dent potential proposed in [19] 



V(x, t) — aix* + a2(l — rt)x 2 



(4.13) 



For t < 1/r, the potential has a single minimum, for 
t > 1/r, it evolves into a double- well. We consider the 
time interval < t < T — 2/r which places the transition 
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(1)) 




(c) 




Fig. 4.1: Numerical determination of the asymptotic work distribution for the evolving double well (|4.13l) with a\ = 1/2, 
a 2 = 6, T — 1 and f3 = 1: (a) optimal trajectories x(t; W), colours code values of the work W, (b) Lagrange parameter 

iq, (c) determinant det A of the fluctuation operator A, and (d) quadratic form { V^' | ^4 1 1 V^) , all as function of W. 



at t — T/2. In contrast to the previous example, neither 
the work distribution, nor its asymptotics can be deter- 
mined using solely analytical techniques. We will therefore 
generate the work distribution from simulations and solve 
the equations fixing the asymptotics numerically. 
To begin with, we have from (|2.3[) and (|2.7[) 

Z = Jdx exp [ - /3(aiir 4 + a 2 x 2 )} (4.14) 
M = exp [(6a 2 x 2 + a 2 )T - ^a 2 rT 2 ] . (4.15) 

The ELE (12TT4]) reads 

x =48ai^ 5 + 32aia 2 (l - rt)x 3 

+ [4o|(l - rt) 2 + 2a 2 r{\ - iq)]x , (4.16) 



and its boundary conditions (|2.15[) are of the form 



x 
x T 



4aiXg 



2a 2 x , 



-4aix T - 2a 2 (l 



rT)x T 



The constraint (|2.4p is 



W = a 2 r I dtx(t;qy 
'o 



(4.17) 



(4.18) 



These equations can be solved numerically using a stan- 
dard relaxation algorithm. The resulting optimal trajecto- 
ries x(t; W) and the corresponding Lagrange parameters 
iq(W) are shown in Fig. I4.1al and Fig. I4.1b[ respectively. 
Due to the mirror symmetry of the potential, there are for 
each value of W two optimal trajectories ±x(t; W), from 
which we only display the positive one. The operator A 
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1 simulation 
- asymptotes 
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•4 
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-20 



-15 

w 



(a) 



10 
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O simulation (10 
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14 
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10 



x 10 



x 

o 



a, 4 



-40 



i 



1 



[ | simulation (10 ) 
I | simulation (10 8 ) 
asymptote 



l 




-35 



-30 



-25 

w 



-20 



-15 



-10 



(c) 



O simulation (10 
X simulation (10 s 
— asymptote 




Fig. 4.2: Work distribution for the evolving double-well (I4.13j) for a\ — 1/2, a.2 — 6, T — 1 and (3 = 1. The histogram 
and the symbols show results from simulation of the Langevin dynamics (|2.1[) . the lines give the asymptotic forms 
(15.301) (full) with and (|2.17p (dashed) without the pre-exponential factor. Subfigures (a) and (b) show a linear and 
logarithmic plot respectively of the work distribution itself, subfigures (c) and (d) display the distribution weighted 
with the factor e~^ w as appearing, e.g., in the Jarzynski equality In (b) and (d) circles and crosses represent 

histograms based on 10 5 and 10 s work values respectively, in (c) results from 10 8 trajectories are shown in light blue. 



from (|2.20p acquires the form 
d 2 

A = — — + 240a?x 4 + 96aia 2 (l - rt)x 2 
at 2, 

+ 2ct 2 r(l - iq) + 4al(l - rt) 2 (4.19) 
with the boundary conditions (|2.21[) 

[I2 ai x 2 + 2a 2 ] tp n (0) - <p n (fl) = , 
[I2ai4 + 2o 2 (l - rT)]<p n (T) + <p n (T) = . (4.20) 
To obtain detA from p.3p . we determine according to 

(32) 

F(0) = [U ai x 2 T + 2a 2 (l - rT)] X o(T) + ±o(T) (4.21) 



by solving numerically the initial value problem (13.81) 

X Q (t) =[240a 2 x 4 + 96aia 2 (l - rt)x 2 

+ 2a 2 r(l - iq) + 4a 2 2 (l - rt) 2 ] X Q(t) = 0, 
*o(0) = 1, Xo(T) = 12ar4 + 2a 2 . (4.22) 

Note that this has to be done for each value of W sep- 
arately by using the appropriate results for x(t; W) and 
q(W). The result for det A as a function of W is depicted 
in Fig. I4.1cl 

The last ingredient for the pre-exponential factor is 
(y'|j4 _1 |y) from p.2j) . To determine it, we need to solve 
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the boundary value problem (|3.1[) . (|2.21l) 

ij)(t) - [240a?x 4 + 96aia 2 (l - rt)x 2 

+ 2a 2 r(l - iq) + 4a|(l - rt) 2 ] ip(t) + 2a 2 rx 
^(0) = (l2a 1 x 2 +2a 2 )i/j{0) 

rj>(Q) = - [l2a lX 2 T + 2a 2 (l - rT)} ip{T) (4.23) 
for each x(t; W) and q{W) and use the result in (|3.2I) 

= -2«2r / dt x(t)ip(t) . (4.24) 
Jo 

The values for (V"'|j4 — -"^IV") obtained in this way are shown 
in Fig. l4~Tdl 

Plugging the numerical results for Af, Zq, x, iq, detA 
and into (12.291) and adding an additional fac- 

tor 2 to account for the two equivalent solutions ±x(t; W) 
for each value of W, we obtain the final result for the 
asymptotic form of the work distribution. 

To investigate the accuracy of this result, we employed 
the Heun scheme to simulate the Langevin equation (|2.1j) . 
In Fig. I4.2al we show the resulting histogram of 10 8 work 
values and the asymptotic behaviour determined above. 
The dashed lines represent the incomplete asymptotic form 
(|2.1 7|) without pre-exponential factor, whereas the full line 
shows the complete asymptotics. In the former case an 
overall constant factor has to be adjusted, in the latter 
no free parameters remain. If the region of work values 
—30 < W < —20 accessible from the simulation using 10 8 
trajectories is utilized for the fit in the incomplete asymp- 
totics, the two asymptotic expressions almost coincide in 
the tail of the distribution. Away from the asymptotic 
regime, however, they differ markedly from each other. If, 
therefore, less data would be available, the fitted incom- 
plete asymptotics could badly fail to reproduce the true 
asymptotic behaviour, see also Fig. I5.1bl The parameter- 
free complete asymptotics is clearly advantageous. 

Furthermore, there is a broad range of excellent agree- 
ment between histogram and complete asymptotics. This 
becomes in particular apparent when examining the 
weighted work distributions P(W) exp(— f3W) shown in 
Fig. l4.2cl and Fig. l4.2dl The average (exp(— (3W)) appear- 
ing in the Jarzynski equality (jl.ll) could already be accu- 
rately determined without the histogram at all by using 
nothing more than the complete asymptotics of P(W). 



5 The breathing parabola 

A particularly interesting class of examples is provided by 
harmonic oscillators with time dependent frequency [201 

V(x,t) = ^x 2 . (5.1) 

Except for some special choices of k(t), the full pdf of 
work is not known analytically. For our purpose the case 



of a monotonously decreasing function k(t) is most ap- 
propriate. Then W < and we aim at determining the 
asymptotic form of P(W) for W — > -co. 
The ELE fiTTty is given by 

x+ ((l~iq)k-k 2 )x = (5.2) 

whereas the boundary conditions (|2.15[) acquire the form 

xq = koXo and xt = — kx xt ■ (5-3) 

These equations constitute themselves a Sturm-Liouville 
eigenvalue problem. Consequently, there are infinitely 
many values q^°\ q^, ■■■ for q. Due to the mirror sym- 
metry of the potential, each value q^ again admits two 
non-trivial solutions ±x' n )(-). 

The somewhat unusual feature of this situation is that 
the different solutions q^ are not related to the value W 
of the constraint. Also, the functional form of x^ n '(-) is 
independent of W. The connection with the work value 
W under consideration is brought about exclusively by 
the pref actor of $W(<) which in view of (|2.4[) and (|5.1I) 
must be 

For the saddle-point approximation in (|2.8[) this means 
that for any value of W there are infinitely many station- 
ary points of P[x(-)]. Moreover, from f|2 . 20[) and (12.21[) we 
find 

A = -^-(l-iq)k + k 2 , (5.4) 
fco Xq — ±o = , kf Xt ~\~ Xt = . (5-5) 

Comparing this with (I5.2j) . (|5 .3[) it is seen that the fluctu- 
ations around the optimal path x(-) are governed by the 
same operator as the optimal path itself, as usual for a 
quadratic action. Together with the homogeneous bound- 
ary conditions (I2.15[) . this implies that every saddle-point 
{q( n \ ±x(™)(-)} has a zero mode, namely the optimal path 
jWf) itself. From Courant's nodal theorem we know that 
has n nodes. But if the Hessian of the saddle- 
point {<f n \ ±x(™)(-)} has the eigenfunction x^ n '{-) with 
zero eigenvalue and n nodes, it consequently must have 
(rt — 1) eigenfunctions with negative eigenvalues. There- 
fore, only the solutions {q^°\ ix^O)} of the ELE corre- 
spond to maxima of P[x(-)], all other solutions are saddle- 
points with unstable directions. These solutions are irrel- 
evant for the asymptotics of P(W), and it is therefore 
sufficient to determine the solutions to (I5.2p . (|5.3p with 
the smallest value g^°^of q. To lighten the notation, we 
will denote the corresponding solutions in the following 
simply by {q, ±x(-)}- 

The general expression (|2.18[) for S greatly simplifies 
for a parabolic potential. First xV — 2V, and the second 
and the third term in (|2.18p cancel. Similarly, xV" = V 
and the forth term vanishes. Finally, xV' = 2V and the 
last term becomes proportional to W. We hence find the 
compact expression 

S=^\W\. (5.6) 
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A complication also arises in the determination of the 
pre-exponential factor. Here we cannot use (|2 . 29[) because 
the zero-mode makes det A — and hence A -1 becomes 
singular. Nevertheless, we may proceed as in section [5] up 
to eq. (|2.23[) which now reads 



I 



n 



dc„ 



dr 



exp 



I 
4 



n>l 



^ / j A n c n 2 / j C ™ a 



(5.7) 



Integrating over c n with n > 1 yields 
/Sir 1 f dco 



/? VdekA 1 J 
f dr / fir 2 s--^ d 2 ij3r 



where we have used the usual notation 
det A' := A„ 



(5.8) 



(5.9) 



for a determinant omitting the zero mode. The remaining 
integrals over r and Co are Gaussian and give 



I : 



V2 



^d 2 det A 1 



(5.10) 



so that we end up with 



P{W) = 2 J -J- (i + (1/ f})) , (5.11) 

Zo y/dl det A' 

where the leading factor of 2 again accounts for the two 
equipollent saddle-points {q,±x(-)}. Comparing this re- 
sult with (|2.29[) . we realize that in the presence of a zero- 
mode we have to replace det^4 by det A' d%. 
In view of (|2.28[) . this is quite intuitive: With Aq tending to 

zero, (U'|A _ \V') becomes more, and more dominated by 

do/Ao and cancelling Ao between det A and (U'|A _1 |U') 
leaves us with (|5.1ip . 

We may finally express do and det A' in terms of x. The 
former is calculated from its definition (|2.24[) and (po — 



d 



i fx- 2W 

„ dtwy= /„ d< H te = H' (5 ' 12) 



The determination of det A' may be accomplished by im- 
plementing a slight variation of the method described in 
section [3] [7| . As shown in appendix [Cj eq. (|3.3[) has to be 
replaced by 

detA ' - # (°) (5.13) 



where 



det A T , 



F(0) 



F rof (0) 

d\\xo\\ 2 
Xo(T) ■ 



(5.14) 



Using d — 1 and xo — x/x(0), we find from (|3.10l) 
det^ = 2Ml. 



X X T 

Therefore, (|5.11l) may be written as 



(5.15) 



P W = ^^ |ff| ( 1 + W))' ^ 
Af = exp(~^ di*(t)) (5.17) 

(5.18) 



where 



and 



are easily calculated. 

We hence find for parabolas with time dependent fre- 
quency the universal asymptotic form 



P --pc 2 \w\ 



(5.19) 



with only the constants C\ and C2 depending on the spe- 
cial choice for kit). Note also that in this case the solution 
of the ELE is sufficient to get the full asymptotics includ- 
ing the prefactor. 

As a simple example we first discuss the case 

t HMt- <«») 

Here, P(W) may again be calculated exactly. The parti- 
cle gains energy only at t = r, where its position is still 
distributed according to po(x). With A — k — kr > 0, we 
hence find 

P(W)= f dx T Po (x T ) 5(W + -x 2 T ) 



exp 



(5.21) 



y nA\W\ 
The ELE (J5HJ is of the form 

x - (1 - iq)A5(t - t)x - k 2 x = (5.22) 
and has the solution 




m\ e ko(t-r) for0 <<<r 

-A- (5.23) 

mi e -k T (t- T ) i0YT<t<T 



where 



(!-»?) = - 



A 



(5.24) 



Using this result in (|5. 16|) and (|5.6p . we find back (|5.21|) . 
A somewhat more general example is given by 



k(t) = 



1 

T+t 



(5.25) 
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Fig. 5.1: Work distribution for a breathing parabola (|5.1I) with protocol (|5.25[) for T — 2 and (3 = 1. The histogram 
and the symbols show results from simulation of the Langevin dynamics (|2.1I) . the lines give the asymptotic forms 
(|5.30D (full) with and (|2.17p (dashed) without the pre-exponential factor. Subfigures (a) and (b) show a linear and 
logarithmic plot respectively of the work distribution itself, subhgures (c) and (d) display the distribution weighted 
with the factor e~^ w as appearing, e.g., in the Jarzynski equality (jl.ll) . Circles and squares represent histograms 
based on 10 5 and 10 8 work values respectively, in (c) results from 10 s trajectories are shown in light blue. 



Here P(W) is not known exactly, nevertheless, some an- Here, with v = 2/xln(l + T) 
alytical progress can be made in the determination of its 
asymptotics. To begin with, we have 



M=VT+T 



and 



Zn = 



1, , 1, 
a ) sin v — cos v + 1 + via H 



(5.26) 



The ELE (15.21) can be solved analytically with the result 

W\ 



and [i = y/iq — 9/4 is the smallest root of 
(4^2 — 3) sin — — 8fi cos — = . 



>0, 
(5.28) 

(5.29) 



x(t) 



± 



Vi + t 



This equation has to be solved numerically. The solution 
for fj, yields the value of iq, determines the prefactor of x 



(2/i cos( A1 ln(l + t)) + sin(/i ln(l +*))). (5.27) J^ 3 ' and thcrefore fixes the com P lete asymptotics via 
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For the special case T = 2 we find /i 
iq = 3.654 which results into 



1.184 implying A Calculation of J 



P{W) ~ 1.021 




-1.8270 \W\ 



(5.30) 



In Fig. 15. II we compare this asymptotics with results from 
numerical simulations of the Langevin equation. Here the 
incomplete asymptotics represented by the dashed line 
was fitted to the numerical results from 10 5 realizations. It 
is clearly seen in the logarithmic plots (b) and (d) that this 
procedure does not yield reliable results for the far tail of 
the distributions. The full asymptotics including the pre- 
exponential factor is again clearly superior and describes 
the true distribution up to rather large values of W. 



6 Conclusion 

In the present paper, we have shown how the complete 
asymptotic behaviour of work distributions in driven 
Langevin systems may be determined. The calculation of 
the pre-exponential factor was accomplished by a method 
building on the spectral (^-function of the operator de- 
scribing the quadratic fluctuations around the optimal 
trajectory. We have shown that the inclusion of the pre- 
exponential factor improves the asymptotics significantly 
and simplifies its application due to the absence of free pa- 
rameters. For the examples considered, our results for the 
asymptotics match the outcome of extensive simulations 
perfectly and reach far into the region of work values acces- 
sible to simulations or experiments with moderate sample 
sizes. 

For the class of harmonic systems with time-dependent 
frequency we established the universal form 



P(W) 



(6.1) 



of the asymptotic work distribution with only the two con- 
stants C\ and Ci depending on the detailed time-depen- 
dence of the frequency. We note that this form is at vari- 
ance with the findings of Speck and Seifert |22) who claim 
that the work distribution for Langevin systems with slow 
but finite driving must always be Gaussian. Our general 
result (|6.1[) shows that while being presumably correct 
for the central part of the distribution, their statement 
does not hold for the asymptotics. We also note that the 
same asymptotic behaviour was found recently for a two- 
dimensional Langevin system with linear non-potential 
forces [23] . 

In view of the complex general expression (|2.18p for the 
exponent in the asymptotics, in which all terms depend 
in a rather implicit way on the work value W, a further 
identification of universality classes for asymptotic work 
distributions remains an open challenge. 

We would like to thank Daniel Grieser for clarifying remarks 
on the determinant of Sturm-Liouville operators and Markus 
Niemann for helpful discussions. 



Consider as special case of (|2.ip the Langevin equation 
X = -CX + y/2//3 £(t) (A.l) 

with some real constant c. According to (I2.6P and (|2.7[) the 
propagator of the corresponding Fokker-Planck equation 
is given by 



p(x T ,T\x o ,0) 

x(T) 
p cT/2 



Vx{-) exp ( - £ / di(i + cx) 2 ) . (A.2) 



s(0)=so 

From the normalization condition J dxr p(%t, T\xq, 0) = 1, 
we have 



1=e cT/2 IW±A 
V 47T 



dxo exp 



/3(a + c) 2 



5) 



x(T)=x T 



x J dx T J T>x(-) exp ^ — — J dt (x + exf^j , 

x(Q)=x 

where a > —c denotes some other real constant. By partial 
integration we find 



1 =e 



;T/2 



P(a + c) 



4tt 



dxo / dxr 



x(T)=x T 

x J Vx 

x(0)=x 



(-)exp(- 



(ax — x )x + (cx T + Xt)xt 



r T d 2 
+ / dt x( - --^+c 2 )x 



dt 2 



t/2 0(a + c) 



J 



1 



4?r y/detA 
where the operator A is defined by 



(A.3) 



Aip = —<p + c 2 ip , 

a<p(0)-<p(Q) = 0, ctp(T) + tp{T) = . (A.4) 
With the methods described in section [3] we easily find 

det A = 2(a + c) e cT , (A.5) 
and comparison with (| A.3|) yields 

(A.6) 

Note that this value does not depend on a and c and is 
therefore valid for all situations considered in the present 
paper. 
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B Calculation of det A vci and F re f(0) 

In this appendix we calculate the determinant of the op- 
erator 

d 2 

Aref = -— (B.l) 



dt 2 



with boundary conditions 



a<p(t = 0) + b(p(t = 0) =0 
c(p(t = T) + d<p(t =T) =0 . 



(B.2) 
(B.3) 



The main point for introducing this reference operator is 
that the ratio of determinants of two Sturm-Liouville oper- 
ators has much nicer analytic properties than each deter- 
minant individually |14j . Intuitively, this is related to the 
fact that in the naive interpretation of the determinant as 
product of all eigenvalues each determinant separately is 
clearly infinite, whereas their ratio may remain bounded. 
Our derivation closely follows [14] . where the correspond- 
ing analysis for Dirichlct boundary conditions was given. 
We start with 



where 



det A-cf = e" c "f (0) 



Cref(s) := J2 A rcf% 



(B.4) 



(B.5) 



is the spectral ^-function of the operator A Te { with A rc f,„ 
denoting its eigenvalues [JJ. Let x\(t) be the solution of 



Xx(t) = X X x(t), Xa(0) = -b, x A (0) = a . (B.6) 



Then 



F ref (A) := c X x(T) + d X x(T) 



(B.7) 



has zeros at all eigenvalues A = A rc f jn of A le {. Correspond- 
ingly, dlni^ re f(A)/dA has poles at these eigenvalues with 
residua given by their multiplicities. We may hence write 



Cef(s) 



1 

2?ri 



dAA- s — l nJ F ref (A) 
dA 



(B. 



where the contour c starts at A = oo + ie, goes down 
parallel to the real axis, makes a half-circle around the 
lowest eigenvalue A rc f.i > and continues parallel to the 
real axis to A = oo — ie. 

For the simple case of A TC f we may solve (|B.6|) and 
determine F TC f(X) explicitly: 

F ref (A) = (ad - be) cos VXT + (en + bd\) — . (B.9) 



VA 



From this result we find 
d 



dA 



lnF rof (A) ~ VX 



for large |A|. Provided s > 1/2, we may hence deform c 
to the contour starting at A = — oo + ie, going up to a 
half-circle around A = and running back to — oo — ie. 
Taking into account that the integrand has a branch cut 



along the negative real axis and substituting A 
we find (for s > 1/2) 



CrefO) 



sm7rs 



dxx s — lnF ro f (— x) 
ax 



dx x 

— In 

da; 



(ad~bc) cosh(v / xt) 
sinh yfxt 



-x ± ie, 



(B.ll) 



(B.12) 



(ca—bdx)- 



In order to use this expression in (|B.4|) it has to be analyt- 
ically continued to s — 0. No problems arise at the lower 
limit of integration, x — 0. Using the representation 



-lnF rof (-x) 



T 

+ 
d 1 

+ — In 

dx 



1 

2x 

[ad — be 
2y/x~ 

ca/x 



(B.13) 



(1 



bd 



(1 



we see, however, that the first two terms entail divergences 
at large x for s — > 0. Splitting the integral in (|B.12|) at 
x = 1, these dangerous terms may be integrated explicitly 
and the continuation to s = presents no further prob- 
lems. 

At the end we find 



Cref ( s = 0) = — In ( — (cb — ad — caT) 
1 bd 



implying 



det A r , 



— (cb — ad — caT) 



(B.14) 



(B.15) 



To use this result in (|3.3p . we have finally to determine 
F re f (0). To this end, we need the solution xo(t) of 



AefXo = -Xo = 
which is given by 



Xo(0) = -6, Xo(0)=a (B.16) 



Xo(t) =at-b . 



(B.17) 



Then 



(0) = c X o(T) + d Xo (T) = acT - cb + da (B.18) 



which, of course, also follows from (|B.9|) for A — > 0. Com- 
bining (|B~15|l and (|B~18| we arrive at (|3TT0|) . 



(B .10) C Calculation of det A' 



In this appendix we sketch the calculation of the deter- 
minant of a Sturm-Liouville operator omitting its zero 
mode for the case of Robin boundary conditions, where 
we closely follow [7]. The main idea is to replace F(X) as 
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defined in (13.71) by a function that is zero only at the non- 
zero eigenvalues of A and behaves asymptotically in the 
same way as F(X). Then all calculations may be done as 
before, and we end up with a relation similar to (|3.3[) . 

To get an idea how the replacement of F(X) may look 
like, it is instructive to consider the scalar product be- 
tween the zero-mode xo of A and a general solution x\ 
of (13.81) . Note that then \o fulfils both the boundary con- 
ditions at t = and t = T whereas \\ fulfils for general 
A only the one at t = as specified in (|3.8I) . We find by 
partial integration 

Mxo\x\) = (xo\Ax\) 

= -XoXaIo + XoX\\o + (^XoIxa) 

= -^^(c X x{T) + dxx{T)) (C.l) 



16. G. M. Wang, E. M. Sevick, E. Mittag, D. J. Searles, and 
D. J. Evans, Phys. Rev. Lett. 89, (2002) 050601 

17. R. van Zon and E. G. D. Cohen, Phys. Rev. E67, (2003) 
046102, Phys. Rev. E69, (2004) 056121 

18. E. G. D. Cohen, J. Stat. Mech., (2008) P07014 

19. S. X. Sun, J. Chem. Phys. 118, (2003) 5769 

20. C. Jarzynski, Phys. Rev. E56, (1997) 5018 

21. D. M. Carberry, J. C. Reid, G. M. Wang, E. M. Sevick, D. 
J. Searles, and D. J. Evans, Phys. Rev. Lett. 92, (2004) 
140601 

22. T. Speck and U. Seifert, Phys. Rev. E70, (2004) 066112 

23. C. Kwon, J. D. Noh, and H. Park, arXiv : 1102 . 2973 



implying 



F{\) = cxa(T) + d X x{T) = — A( X o|xa) . (C.2) 

Xo(T) 



Hence 



F { A) = ±-±F(\) = (l-\) d{xolxx) 



A 



Xo(T) ' 



(C3) 



vanishes at all eigenvalues A„ > (since F does), remains 
non-zero for A = A = (cf. IC.2[) and behaves asymp- 
totically for large A exactly as F. Similarly to the case 
without zero modes, one hence finds 



detA' 



F(0) 



detA rcf F ref (0) 
which coincides with (|5 . 13[) . 



(C4) 



References 

1. C. Jarzynski, Phys. Rev. Lett. 78, (1997) 2690 

2. D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Phys. 
Rev. Lett. 71, (1993) 2401 

3. G. Gallavotti and E. G. D. Cohen, Phys. Rev. Lett. 74, 
(1995) 2694 

4. U. Seifert, Eur. Phys. J. B64, (2008) 423 

5. M. Esposito and C. Van den Broeck, Phys. Rev. Lett. 104, 
(2010) 090601, Phys. Rev. E82, (2010) 011143 and 011144 

6. D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco 
and C. Bustamante, Nature 437, (2005) 231 

7. K. Kirsten and A. J. McKane, Ann. Phys. (N.Y.) 308, 
(2003) 502 

8. K. Sekimoto Prog. Theor. Phys. Supp. 130, (1998) 17 

9. M. Chaichian and A. Demichev, Path integrals in Physics 
(IOP Publishing, London, 2001) 

10. H. Touchette, Phys. Rep. 478, (2009) 1 

11. I. M. Lifshitz, Sov. Phys. Usp. 7, (1965) 549 

12. B. I. Halperin and M. Lax, Phys. Rev. 148, (1966) 722 

13. A. Engel, Phys. Rev. E80, (2009) 021120 

14. K. Kirsten and P. Loya, Am. J. Phys. 76, (2008) 60 

15. O. Mazonka, C. Jarzynski, arXiv: cond-mat/9912121 



